Final Project

Pandemic Simulator

Team members: Matt Mitchell, Nick Vino, and Logan Jolicoeur

Section I

Introduction

In this notebook, we'll develop a model of a pandemic as it spreads in a susceptible population, and explore the effectiveness of interventions such as vaccination.

According to 2022 census records, the city of Banner Elk in North Carolina has a population of approximately 1,000 people. Banner Elk will be the center of our attention as we simulate the spread of an infection in the city. We introduce well-known models of infectious disease, the SI model, the Kermack-McKendrick (SIR) model, SIRS model, SIRV model and use them to explain the progression of the disease over the course of 2,000 days.

Section II

Implementation

The following models were produced with the Epidemics on Networks (EoN) package for python

Citations

Journal of Open Source Software publication

EoN documentation

The following packages are necessary for this notebook:

import EoN
import networkx as nx
from collections import defaultdict
import matplotlib.pyplot as plt
import numpy as np
from PIL import Image
from collections import Counter
from IPython.display import HTML, display, Markdown, Latex

SIS Model

This model consists of nodes (individual) and edges (links between nodes).

Each node has an indepedent stage susceptible and infected. Any susceptible node(person) if infected, will again become susceptible and may be infected again at some point. Infection can take place at any time in this model.

Initially, we consider a model with no recovery rate. The susceptible compartment will flow entirely into the infected compartment over time.

This model is the most basic of all models in this document.

We can set parameters of the model as follows:

N = 1000            # Set the population
rho = 0.01          # Intial infected (10)
kave = 5            # Average degree
gamma = 0           # Recovery rate
tau = 0.002         # Transmision rate

S0 = (1-rho)*N
I0 = rho*N
SI0 = (1-rho)*kave*rho*N
SS0 = (1-rho)*kave*(1-rho)*N
t, S, I = EoN.SIS_homogeneous_pairwise(S0, I0, SI0, SS0, kave, tau, gamma, tmax=2000)

plt.plot(t, S, label = 'S')
plt.plot(t, I, label = 'I')
plt.xlabel('$t$')
plt.ylabel('Population')
plt.legend()
plt.show()

SIR Model

The SIR model is similar to the SIS model, but in this model after susceptible become infected they become R(t) or recovered. Where they cannot become infected again.

The susceptible equation is as follows:

dS/dt = -b * s(t) * I(t)

The recovered equation is:

dr/dt = k i(t)

The infected equation is:

ds/dt + di/dt + dr/dt = 0

These equations are utilized to develop our models below, we some variation.

N=1000         #population
gamma = 0.0045  #recovery rate
tau = 0.007     #transmission rate
kave = 5        #people we know
rho = 0.01      #initial infected
phiS0 = 1-rho
def psi(x):
    return (1-rho)* np.exp(-kave*(1-x))
def psiPrime(x):
    return (1-rho)*kave*np.exp(-kave*(1-x))

t, S, I, R = EoN.EBCM(N, psi, psiPrime, tau, gamma, phiS0, tmax = 2000)


plt.plot(t, S, label = 'S')
plt.plot(t, I, label = 'I')
plt.plot(t, R, label = 'R')
plt.xlabel('$t$')
plt.ylabel('Population')
plt.legend()
plt.show()

As seen above we follow a typical SIR Model.

SIRS Model

The SIRS Model is a model based on the SIR, however, the recovered will eventually become S(t) again. This can add variation in any simulator, as persons can become infected again, after an infection has occured in a population region.

N = 1000
kave = 5    #expected number of friends/family
G = nx.fast_gnp_random_graph(N, kave/(N-1))

H = nx.DiGraph()
H.add_edge('I', 'R', rate = 0.0025)   #Recovery rate
H.add_edge('R', 'S', rate = 0.001)   #Recovered -> Susceptible

J = nx.DiGraph()  
J.add_edge(('I', 'S'), ('I', 'I'), rate = 0.0065)  #IS->II

IC = defaultdict(lambda: 'S')
for node in range(10):
    IC[node] = 'I'

return_statuses = ('S', 'I', 'R')

t, S, I, R = EoN.Gillespie_simple_contagion(G, H, J, IC, return_statuses, tmax = 2000)

plt.plot(t, S, label = 'S')
plt.plot(t, I, label = 'I')
plt.plot(t, R, label = 'R')
plt.ylabel('Population')
plt.xlabel('$t$')
plt.legend()

<matplotlib.legend.Legend at 0x1e316070>

SIRV Model

G = nx.grid_2d_graph(32,32) #each node is (u,v) where 0<=u,v<=32
#we'll initially infect those near the middle
initial_infections = [(u,v) for (u,v) in G if 14<u<18 and 14<v<18]
H = nx.DiGraph() #the spontaneous transitions
H.add_edge('Sus', 'Vac', rate = 0.0003)
H.add_edge('Inf', 'Rec', rate = 0.0025)
J = nx.DiGraph() #the induced transitions
J.add_edge(('Inf', 'Sus'), ('Inf', 'Inf'), rate = 0.01)

#calculate R for the graphs
def calc_R(v):
    check_stat = []
    R = 0
    for time, w, *h in v:
        check_stat.append([w])
     
    c = Counter(i[0] for i in check_stat)
    for x, y in c.items():
        if(y > 0 and len(c.items()) > 0):
            R += y 
            
    return R / len(c.items())
    
IC = defaultdict(lambda:'Sus')
for node in initial_infections:
    IC[node] = 'Inf'
    
return_statuses = ['Sus', 'Inf', 'Rec', 'Vac']

color_dict = {'Sus': '#001aff','Inf':'#ff2000', 'Rec':'#008000','Vac': '#fffb00'}
pos = {node:node for node in G}
tex = False
sim_kwargs = {'color_dict':color_dict, 'pos':pos, 'tex':tex}

sim = EoN.Gillespie_simple_contagion(G, H, J, IC, return_statuses, tmax=2000, return_full_data=True, sim_kwargs=sim_kwargs)

display(Markdown(str(calc_R(sim.transmissions()))))

times, D = sim.summary()

newD = {'Sus+Vac':D['Sus']+D['Vac'], 'Inf+Rec' : D['Inf'] + D['Rec']}

new_timeseries = (times, newD)
sim.add_timeseries(new_timeseries, label = 'Simulation', color_dict={'Sus+Vac':'#E69A00', 'Inf+Rec':'#CD9AB3'})

sim.display(time=1, node_size = 4, ts_plots=[['Inf'], ['Sus+Vac', 'Inf+Rec']])

ani=sim.animate(ts_plots=[['Inf'], ['Sus+Vac', 'Inf+Rec']], node_size = 4)
ani.save('SIRV_animate.gif', writer="pillow", dpi=100)

1.46875

G = nx.grid_2d_graph(32,32) #each node is (u,v) where 0<=u,v<=32
#we'll initially infect those near the middle
initial_infections = [(u,v) for (u,v) in G if 14<u<18 and 14<v<18]
H = nx.DiGraph() #the spontaneous transitions
H.add_edge('Sus', 'Vac', rate = 0.0005)
H.add_edge('Inf', 'Rec', rate = 0.0025)
J = nx.DiGraph() #the induced transitions
J.add_edge(('Inf', 'Sus'), ('Inf', 'Inf'), rate = 0.01)

#calculate R for the graphs
def calc_R(v):
    check_stat = []
    R = 0
    for time, w, *h in v:
        check_stat.append([w])
     
    c = Counter(i[0] for i in check_stat)
    for x, y in c.items():
        if(y > 0 and len(c.items()) > 0):
            R += y 
            
    return R / len(c.items())

persons = {}

IC = defaultdict(lambda:'Sus')
for node in initial_infections:
    IC[node] = 'Inf'
    
return_statuses = ['Sus', 'Inf', 'Rec', 'Vac']

color_dict = {'Sus': '#001aff','Inf':'#ff2000', 'Rec':'#008000','Vac': '#fffb00'}
pos = {node:node for node in G}
tex = False
sim_kwargs = {'color_dict':color_dict, 'pos':pos, 'tex':tex}

sim = EoN.Gillespie_simple_contagion(G, H, J, IC, return_statuses, tmax=2000, return_full_data=True, sim_kwargs=sim_kwargs)

display(Markdown(str(calc_R(sim.transmissions()))))

times, D = sim.summary()

newD = {'Sus+Vac':D['Sus']+D['Vac'], 'Inf+Rec' : D['Inf'] + D['Rec']}

new_timeseries = (times, newD)
sim.add_timeseries(new_timeseries, label = 'Simulation', color_dict={'Sus+Vac':'#E69A00', 'Inf+Rec':'#CD9AB3'})

sim.display(time=1, node_size = 4, ts_plots=[['Inf'], ['Sus+Vac', 'Inf+Rec']])

ani2=sim.animate(ts_plots=[['Inf'], ['Sus+Vac', 'Inf+Rec']], node_size = 4)
ani2.save('SIRV_animate.gif', writer="pillow", dpi=100)

1.4351585014409223

HTML(ani.to_jshtml())

Once Loop Reflect

HTML(ani2.to_jshtml())

Once Loop Reflect

Masks

In this scenario, we implement masks. We will be using the SIR model to demonstrate how effective masks can be in curtailing the pandemic.

Based on mayoclinic.org report (https://www.mayoclinic.org/diseases-conditions/coronavirus/in-depth/coronavirus-mask/art-20485449).

The indicate that 83% of masks assist with curtailing COVID-19. This is not a representation of how COVID-19 affects the population, but it is a starting point for how effective masks can be in a "general" pandemic scenario. Without masks, we start an infection rate of 0.0065. Taking 83% of 0.0065, we get 0.001105. As such this is the representation of how masks can affect a pandemic, if the population were to wear them at time(1).

N = 1000
kave = 5    #expected number of friends/family
G = nx.fast_gnp_random_graph(N, kave/(N-1))

H = nx.DiGraph()
H.add_edge('I', 'R', rate = 0.0025)   #Recovery rate
H.add_edge('R', 'S', rate = 0.001)   #Recovered -> Susceptible

J = nx.DiGraph()  
J.add_edge(('I', 'S'), ('I', 'I'), rate = 0.001105)  #IS->II 83% of 0.0065 (this is based on the stat that 83% of masks curtail people getting infect with COVID)

IC = defaultdict(lambda: 'S')
for node in range(10):
    IC[node] = 'I'

return_statuses = ('S', 'I', 'R')

t, S, I, R = EoN.Gillespie_simple_contagion(G, H, J, IC, return_statuses, tmax = 2000)

plt.plot(t, S, label = 'S')
plt.plot(t, I, label = 'I')
plt.plot(t, R, label = 'R')
plt.ylabel('Population')
plt.xlabel('$t$')
plt.legend()

<matplotlib.legend.Legend at 0x175d0730>

Without Masks

This scenario shows what happens if we do not have masks and the population does not wear them at all. The rate will be the full 0.0065 rate. Show casing that there is nothing stopping the infection from spreading betweens friends and family, except the chance of getting the virus.

N = 1000
kave = 5    #expected number of friends/family
G = nx.fast_gnp_random_graph(N, kave/(N-1))

H = nx.DiGraph()
H.add_edge('I', 'R', rate = 0.0025)   #Recovery rate
H.add_edge('R', 'S', rate = 0.001)   #Recovered -> Susceptible

J = nx.DiGraph()  
J.add_edge(('I', 'S'), ('I', 'I'), rate = 0.0065)  #IS->II, half the rate of  demonstrating an half in infection occuring if masks are worn by the population

IC = defaultdict(lambda: 'S')
for node in range(10):
    IC[node] = 'I'

return_statuses = ('S', 'I', 'R')

t, S, I, R = EoN.Gillespie_simple_contagion(G, H, J, IC, return_statuses, tmax = 2000)

plt.plot(t, S, label = 'S')
plt.plot(t, I, label = 'I')
plt.plot(t, R, label = 'R')
plt.ylabel('Population')
plt.xlabel('$t$')
plt.legend()

<matplotlib.legend.Legend at 0x160c9810>

As you can tell above. We have a significant increase of people that acquired the infection, than those that didn't. With a mask stopping 83% of an infection actually occuring for S(t) people, we had a significant decrease. Based on our evidence, masks do something and not anything at all. However, different types of masks and their grade is not included into these scenarios, so do keep that in mind.

Section III

Conclusion

To conclude, the different models above assist in a better understanding of how diseases can spread, and with differing variations we can better understand how lockdowns, social gatherings, quarantines, or masks can either improve or unimprove a pandemic. It is important to note that these models do not take into account all variables that can occur during a pandemic, because humans and viruses can be unpredictable sometimes. Which leaves much to understand.